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Abstract. We present a new particle code for modelling the evolution of galaxies. The code is based on a multi-phase de- 
scription for the interstellar medium (ISM). We included star formation (SF), stellar feedback by massive stars and planetary 
nebulae, phase transitions and interactions between gas clouds and ambient diffuse gas, namely condensation, evaporation, drag 
and energy dissipation. The latter is realised by radiative cooling and inelastic cloud-cloud collisions. We present new schemes 
for SF and stellar feedback. They include a consistent calculation of the star formation efficiency (SFE) based on ISM properties 
as well as a detailed redistribution of the feedback energy into the different ISM phases. 

As a first test example we show a model of the evolution of a present day Milky-Way-type galaxy. Though the model exhibits 
a quasi-stationary behaviour in global properties like mass fractions or surface densities, the evolution of the ISM is locally 
strongly variable depending on the local SF and stellar feedback. We start only with two distinct phases, but a three-phase ISM 
is formed soon consisting of cold molecular clouds, a warm gas disk and a hot gaseous halo. Hot gas is also found in bubbles 
in the disk accompanied by type II supemovae explosions. The volume filling factor of the hot gas in the disk is ~ 35%. The 
mass spectrum of the clouds follows a power-law with an index of a a -2. The star formation rate (SFR) is ~ 1.6Mq yr"' on 
average decreasing slowly with time due to gas consumption. In order to maintain a constant SFR gas replenishment, e.g. by 
infall, of the order 1 Mq yr"' is required. Our model is in fair agreement with Kennicutt's (1998) SF law including the cut-off 
at ~ 10 Mq pc"-. Models with a constant SFE, i.e. no feedback on the SF, fail to reproduce Kennicutt's law. 
We have performed a parameter study varying the particle resolution, feedback energy, cloud radius, SF time scale and metal- 
licity. In most these cases the evolution of the model galaxy was not significantly different to our reference model. Increasing 
the feedback energy by a factor of 4 - 5 lowers the SF rate by ~ 0.5 Mq yr"' and decreasing the metallicity by a factor of ~ 100 
increases the mass fraction of the hot gas from about 10% to 30%. 
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1. Introduction 

The treatment of the instellar medium (ISM) is a crucial part 
in modelling galaxies due to the importance of dissipation, star 
formation (SF) and feedbac k. Single -phase TREE-SPH mod- 
els su ffer from over-cooling JPearce et al. 1999: Thacker et al.l 
I2OO0I) because the multi-phase structure of the ISM is not 
resolved. As a result baryonic matter loses angular momen- 
tum to the dark matter (DM) halos and the disks formed in 
these simulations are too concentrated. Also, SF is overly effi- 
cient. On the other hand, a multi-phase ISM is a pre-requisite 
for the classical grid-based chemo-dynamical (CD) models 
llBurkert et alJ Il992| iTheis et alJ Il992t ISamland etal1ll997t 
ISamland & G erharcf 2003*) but thev lack the geometrical flexi- 
bility of a (3d-)particle code. 

An approach to include a multi-phase t reatment of the ISM 
in a particle code has been presented by ISemelin & Combes! 



(l2002h . who combined Smoothed Particles Hydrodynamics 



(SPH) with a sticky 
'Combes & Gerin 'l985"; 



particle method 
iNoguchil Il988l IPalous et all 



Hydrody i 
dCarlberd 



1984| 
1993t 

multi-phase ISM code of 
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Theis & Hensler 1993). In the 

Isemelin & Combesl l 2002h sticky particles are used to describe 
a cold, dense cloudy medium while the diffuse ISM is modelled 
by SPH. Stars form from cold gas and feedback processes in- 
clude heating of cold gas by supernovae and stellar mass loss 
(including metal enrichment). Cooling leads to a phase transi- 
tion from warm to cold gas. The main advantage of this ap- 
proach is the dynamical treatment of the cloud component: 
clouds are assumed to move on ballistic orbits dissipating ki- 
netic energy in collisions. 

The main purpose of this paper is to present a new ex- 
tended TREE-SPH code with a multi-phase descripti on of the 
ISM a dopted from CD models. Similar to Semelin & CombesI 
(l2002h . this is achieved by a combination of a sticky particle 
scheme and SPH. However, we use a different sticky parti- 
cle scheme based on an individual treatment of single clouds. 
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Fig. 1. The network of processes connecting the different 
gaseous and stellar components. Clouds are the sites of SF. 
Stars return mass and inject energy to the multi-phase ISM. 
The two gas phases exchange mass and momentum by means 
of condensation or evaporation and drag. Dissipation of energy 
is due to radiative cooling and cloud-cloud collisions. 

Advantages of this scheme are the build-up of a cloud mass 
spectrum, the individual handling of cloud collisions and the 
localised treatment of phase transitions between clouds and the 
ambient diffu se gas. Furthermore, we apply a refined SF recipe 
suggested bv lElmegreen & Efremo vl ill997l hereafter EE97). 
This scheme considers star formation within individual clouds 
including local feedback processes. As a result the star forma- 
tion efficiency is locally variable depending on the cloud mass 
and the ambient ISM pressure. 

A detailed description of the model including all the pro- 
cesses is given in Sect. |2l Details on the numerical treatment 
can be found in Sect. |3] The code is applied to an isolated 
Milky-Way-type galaxy (Sect.|3. The results (including a short 
parameter study) are discussed in Sect. [S] Conclusions and fu- 
ture applications of this code are presented in Sect. |6l 

2. The Model 

We use particles to describe all galactic components. Each par- 
ticle is assigned a special type representing either stars, clouds, 
(diffuse) gas or DM. The different components are connected 
by a network of processes including SF, stellar feedback, cool- 
ing, cloud-cloud collisions, condensation, evaporation and drag 
as shown schematically in Fig.[n 

2.1. Stars 

Each star particle represents a cluster of stars formed at the 
same time. The initial ste llar mass distr i bution is given by the 
multi-power law IMF of Kroupaet al. (1993) with the lower 
and upper mass limits of 0. 1 Mq to 100 Mq respectively. 

The stellar life times (depending o n stellar mass and 
metallicity) are approximated according to Raiteri et al.Vl996) 
to give reasonable fits to the Padoya evolutionary trac ks 
jAlongi et alJl993HBressan et alJl993HBertelli et alJl994 in 
the mass range from 0.6 Mq to 120 Mq and for metallicities 
Z = 7 ■ 10-5 - 0.03. 



We assume that stars with masses above 8 Mq end their life 
as type II supernovae (SNII) and stars from 1-8 Mq end as 
planetary nebulae (PNe). The mass fractions of these stars are 
9% and 34% respectively. The life times of the high mass stars 
range from 3 Myr to 39 Myr. The stellar feedback is described 
below. 

2.2. Clouds 

The cloud particles describe a cold, clumpy phase of the ISM, 
i.e. the molecular clouds, which are also the sites of SF. The 
radius r^d of an individual cloud i s given by the following re- 
lation based on observations (.Rivolo & Solomonll988i) : 

. / nicid 

where the constant Add is set to 50. 

The giant molecular clouds represented by cloud particles 
consist of gas at diff'erent densities and temperatures. Since 
we do not resolve this structure we do not define a temper- 
ature or pressure for individual clouds. However, we assume 
that the clouds are influenced by the pressure of the surround- 
ing ISM (see Sect, \2.6i . Radiative cooling within clouds is 
not taken into account. However, we allow for dissipation of 
kinetic energy in cloud-cloud collisions using the dissipation 
scheme of Theis & Hensler (.1993). Two clouds undergo a com- 
plete inelastic collision if they physically collide within the 
next time step. Additionally, their orbital angular momentum 
must be less than the break-up spin of the final cloud and 
the cloud mass must not exceed a critical mass limit (here: 
■^coU.iim = 1 ■ 10^ Mq). In casc of a sticky collision the two 
cloud particles are replaced by one particle keeping the center- 
of-mass data of the former two clouds. 

2.3. Diffuse gas 

In an first approximation the diffuse gas can be described as 
continuous self-gravitating fluid. We model this component us- 
ing an SPH approach. The gas is assumed to obey the ideal 
equation of state with adiabatic index 7 = 5/3: 

^'gas = (r - l)PgasM, (2) 

where Pgas is the pressure, p the density and u the energy per 
unit mass of the gas. The thermal energy evolution is governed 
by cooling and heating processes. For the latter we only take 
the most dominant source, the SNII, into account. Dissipation 
is performed by the metal-dependent cooling function from 
JBohringer & Hensler ( 1989). 

The difiiise gas can have temperatures from a few 10^ K 
up to several 10* K. In the following we will also distinguish 
between warm and hot gas meaning gas at temperatures below 
or above 10'* K respectively. 

2.4. Dark Matter 

In our code DM can be modelled either with particles (live 
halo) or as a static potential. In this work we focus on the 



S. Harfst et al.: Modelling Galaxies with a 3d Multi-Phase ISM 



3 




Pgas > PsF.lim 
mold>MsF,„„, 

'''ia < ti - tcld 

P T //£/// 



m 



//// 




Fig. 2. The diagram shows the SF scheme: (fo) the cloud is in- 
active (no SF); (fi) an embedded star cluster has been formed 
with a local SF efficiency; (fa) the cloud is fragmented by SNII 
energy input. 

late evolution of an isolated galaxy where a static DM poten- 
tial is a reasonable approximation. We use the potential from 
iKuiiken & Dubinskil ( ll995l) which is a flattened, isothermal po- 
tential with a core and a finite extent. 

2.5. Interactions between ttie ISM phases 

The ISM is treated as two dynamically independent gas 
phases. This approac h has been successfully applied in chemo - 
dynamical models dTheis et alJ Il992t ISamland et al* 'l997). 
The two gas phases are hnked by two processe s: condens a- 
tion/ evaporation (C/E) and a drag force. Following lCowie et all 
lll981i) the mass exchange rate by condensation and evapora- 
tion for an individual cloud particle is calculated from' 



meld 



^ gas 1 r eld 



pc 



-2.75 
-2.75 
8.25 ■ 



■ lO^ 

loVo' 



if (To > 1 

if 0.03 < o-o < 1 g s- 
if o-o <0.03 



(3) 



where the parameter ctq is defined as 
"^'^ I I 



^ -'gas \ / rM \ ( "gas V 

1.54- lO^K PC \cm-'l 



(4) 



Tgas and n„.^^ are the local temperature and number density of 
the diffuse gas. The clouds are sub ject to a drag fo rce Fd due 
to ram pressure which is given bv llShuet alJl972h 



Fd = -Co^rr^idPgas VrelVrel, 



(5) 



where pgas is the local gas density and vvei is the velocity of 
the cloud relative to the gas. The constant Co is not well deter- 
mined and it can be in the range 0.1 - 1 . We used Cd = 0. 1 . 

2.6. Star Formation 

The treatment of SF is usually based on the Schmidt Law 
llSchmidJl959l) 



SFR/area = c„ ■ 2" 



(6) 



' We changed the factor for the case ctq > I from -3.75 ■ 10"* to 
-2.75 • 10* in order to avoid a discontinuity at ctq = 1. 



with the gas surface density Eg and n a 1.5. The parameter c„ 
combines a SF time scale and a SFE. A reasonable estimate for 
the time scale is given by the dynamical time. The SFE on the 
other hand cannot be determined from current theories of SF 
and therefore remains a free parameter A constant SFE also 
cannot account for local feedback process es, by this prohibiting 
any self -regulation of the SF. However, iKoppen et alJ d 19951) 
have demonstrated the importance of negative stellar feedback 
which leads to a Schmidt-type strongly self-regulated SF. 

Though the standard approach for SF works well for iso- 



lated galaxies, it has ( 
of interacting galaxies 


lifficulties in reproducin 


g observations 


aviihos et alJl99^ 


1 19931) 


Earned 


1I2OO4I) 



suggested a second shock-induced SF mode. However, this 
means another free parameter in the SF recipe. 

Different to such an approach, we applied a single SF law, 
but with a temporally and spatially variable SFE suggested by 
EE97. They assign each molecular cloud an individual star for- 
mation efficiency which is controlled by the ISM properties 
(mass of the molecular cloud, pressure of ambient diffuse gas) 
and the stellar feedback. According to EE97 the star formation 
is most efficient for massive clouds embedded in a high pres- 
sure ambient medium. A pre-requisite of using the EE97 SFE 
scheme is a multi-phase treatment of the ISM as proposed in 
this paper 

The main concept of our SF implementation is that stars 
form in clouds and clouds are destroyed by stellar feedback 
and, thus, self -regulating the SF (Fig.|5J. Each cloud is assumed 
to form stars only after a time of inactivity Tia (to in Fig. |2}, 
because not all gas in clouds is dense or unstable enough for 
immediate SF. Once the SF process has started, an embedded 
star cluster is formed (ti) instantaneously ^. 

Following |eE97 we use a variable SF efficiency e depend- 
ing on local properties of the ISM. A typical SFE (cloud mass 
of 10^ - 10^ Mq and pressure typical for solar vicinity) is be- 
tween 5% and 10% (see Fig. 4 in EE97). For details of the 
implementation see App.lAl 

The only free parameter in our SF recipe is ria. It can be 
interpreted as a global SF time scale: assuming a SFE of e ~ 
0.1, a total mass Mdd in clouds of a few 10^ M0 and Tia of a 
few 100 Myr , we get a global SFR ~ eMddT;^' of ~ 1 Mq yr\ 

2.7. Stellar Feedback 

Once the SFE has been determined the energy released by SNII 
in each time step can be computed from the chosen IMF and 
stellar life times (Sect. IZTl . We assume esN = 2 ■ lO'^'erg for 
the energy released per SNII event. 

This energy injection disrupts the cloud into smaller frag- 
ments (t2 in Fig. 13 leaving a new stellar particle. The time for 
disruption tdm i- t2 - ti) is determined from the energy input: 
in our model the energy input of the SNII drives an expand- 
ing shell into the cloud. The cloud mass is accumulated in this 
shell. The expansion of the SN shell can be calculated from 



^ We neglect the short time for the SF process itself. This is mo- 
tivated by recent observations suggesting that SF occurs on the local 
crossing time scale for sound waves in a cloud which is shorter than a 
few 10 Myr liElmegreea,200Q) 
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a self-similar solution dBrown et alJll995t iTheis et alJll998l) . 

Assuming a spherical cloud with a 1/r-density profile in agree- 
ment with the mass-radius relation (Eq.0 we get 

/ E \° '' 

rsh(f) = 0.918(-l -f"", (7) 
/ E \°-^' 

Vsh(f) = 0.6891-1 (8) 

where E is the energy injection rate and pi - m^id/r^^^. 

It is assumed that the cloud is disrupted when the shell ra- 
dius reaches the cloud radius. The fragments are then placed 
symmetrically on the shell with velocities corresponding to the 
expansion velocity at fd;,. The energy not consumed by cloud 
disruption, i.e. the energy of SNII exploding after fdis, is re- 
turned as thermal feedback to the hot gas phase. Depending on 
cloud mass and SF efficiency the fragmentation uses about 0.1- 
5% of the total SNII energy of which ~ 20% end up as kinetic 
energy of the fragments. 

The mass returned by SNII is also added to the surrounding 
gas. It is determined from the IMF and stellar life times assum- 
ing^t hat the mass ejected per SNII event MsNn.ej is given by 
fRaiterietalJl996h 

MsNu,ej(m) = 0.77m'0^MQ, (9) 
where m is the mass of a star 

Additionally to the SNII feedback, we also allow for feedback 
by PNe. After each time step, the mass MpN returned to the ISM 
by PNe is calculated by assuming that the PN mass ejected by 
a star of mass m is given bv ( iWeidemaniji2000l) 

Mpn,ej(m) = 0.91m - 0.45 Mq. (10) 

Gradual metal enrichment due to stellar mass release and 
type I SNe are not included so far. This is justified for the kind 
of models presented here (s. Sect. 14. U as metal enrichment has 
little effect as we show later (s. Sect. 15. 3t . 

3. Numerical Treatment 

3.1. Stellar Dynamics 

In our code the gravitational force on a particle can be calcu- 
lated as a combination of the self-gravity of the N-body system 
and an external force. External forces - if applied - usually 
mimic a static dark matter halo. Self-gravity is calculated us- 
ing a tree scheme. We have included two different tree schemes: 
the widely used scheme of.Barnes & Hu t (J_986) ^nd a new tree 
algorithm proposed by IPehnenl J2002l) . An advantage of the 
Dehnen-TREE is that it is momentum-conserving by construc- 
tion. Additionally, an error control is included and the CPU 
scaUng is basically linear with the number of particles. In our 
implementation the Dehnen-TREE is more than an order of 
magnitude faster than the Barnes&Hut-TREE. Hence, we ap- 
plied here mainly the Dehnen-TREE. For easy comparison with 
other N-body solvers we realised gravitational softening by a 
Plummer potential with a softening length es = 0.1 kpc. 



3.2. Gas Dynamics 

We use the SPH formalism to compute the hydrodynamics of 
the diffuse gas and a sticky particle scheme for the clouds. 
Details of our SPH implementation are described in Add. 151 

The dynamics of the clouds are described by the sticky par- 
ticle scheme of Theis & Hensler ( 1993). According to them 
clouds move on ballistic orbits unless they undergo inelastic 
collisions. The treatment of cloud-cloud collisions is described 
in Sect.O 

3.3. Time Integration 

The time integration of the system is done with a leap-frog 
scheme. Individual time steps are used for each particle. The 
time step for an SPH particle p is chosen to satisfy a modified 
form of the Courant condition (for details see App. [EJ. The 
time step of all other particles is based on dynamical criteria 
calculated by 

At,, mini ^, ^P^, (11) 

\vp yap) 

where Op is the acceleration and the softening length of par- 
ticle p. The constant T/ts is of the order of unity (we set it to 
0.4). 

The largest individual time step Atp also defines the system 
time step Afsys which is further limited not to exceed ~ 1 Myr. 
Individual particle ti me steps are c hosen to be a power-of-two 
subdivision of Afsys llAarsethll985h . 

3.4. Mass exchange and new particles 

Several processes included in our model lead to an exchange of 
mass between particles or the creation or removal of particles. 
These processes are introduced in a matter- and momentum- 
conserving manner. 

The mass exchange by condensation and evaporation is cal- 
culated after each system time step Afsys. The transferred mass 
mcidAfsys is added to the cloud. The same mass is subtracted 
from the surrounding SPH particles taking into the account the 
weighting from the kernel function. Additionally, the velocity 
of particles gaining mass (i.e. clouds in case of condensation 
and SPH particles for evaporation) is changed in order to con- 
serve momentum. The mass exchange from stars to SPH parti- 
cles (SNII) and clouds (PNe) is done in a similar way. 

Whenever a SF process occurs, the star forming cloud is 
removed from the simulation. Instead a new star particle of the 
mass emcid is added at the same position with the same ve- 
locity. Furthermore, the clouds formed in expanding shells are 
created by adding four new cloud particles carrying the remain- 
ing mass. Their positions and velocities are chosen according 
to Eq. and (|8} thus conserving mass and momentum. 

3.5. Mass limits for particles 

Our model allows particles to be created, merged and changed 
in mass. A few constraints are used for the sake of numerical 
accuracy and particle numbers. 
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Table 1. Properties of the initial disk galaxy model 



component 


mass' 


radius- 


No. of particles 


disk 


0.34 


18.0 


20000 


stars 


0.27 




16000 


clouds 


0.07 




4000 


bulge 


0.17 


4.0 


10000 


DM halo 


1.94 


87.9 




hot gas 


0.01 


25. 


5000 


total 


2.47 




35000 



' in units of 10" Mq, ^ in kpc 



In order to prevent the formation of too many low mass star 
particles a lower mass limit Msp.iim = 2.5 ■ 10** Mq and a lower 
pressure limit PsF.iim = O.lf for SF are introduced. 

The maximum mass for clouds is given by Mcoiijim 
(Sect. 12. 2^ but a lower mass limit for clouds Mddjini is also set 
to 

Mdd,lim = O.lMsF.lim- (12) 

If the mass of a cloud becomes smaller than Mdd.iim the particle 
is removed from the system and the mass is distributed evenly 
among neighbouring cloud particles within 300 pc. 

SPH particles can change their mass due to condensation 
and evaporation. In order to avoid SPH particles getting too 
massive an upper mass limit is introduced. Whenever an SPH 
particle reaches this limit it is split into four new particles. The 
new equal mass particles have the same center of mass, ve- 
locity, angular momentum and specific energy as the particle 
they replace. A lower mass limit is also used where the mass 
of an SPH particle is distributed to the neighbouring particles. 
The upper and lower mass limits are set to 4 mspH and 0. 1 otsph 
respectively, where mspH is the SPH particle mass at the begin- 
ning of the simulation. 

3.6. Tests 

We have tested our code thoroughly. Energy and angular mo- 
mentum conservation is better than 1 % for a pure A^-body sim- 
ulation following the evolution of a disk galaxy over many ro- 
tational periods. We tested our SPH implementation with the 
standard Evrard-test tevrardl[l98 i) and find that energy con- 
servation is within 3%. In a full model, angular momentum is 
conserved within 3%. 

4. Disk galaxy model 

4.1. Initial Conditions 

A galaxy similar to the present-day Milky Way was set up using 
the ga laxy building package described bv Kuiiken & Dubinski 
lll995[ hereafter KD95). The models of IkDQsI consist of three 
components: the baryonic disk and bulge, and a DM halo. A 
distribution function for each component yields a unique den- 
sity for a given potential. By solving the Boltzmann equation 
for the combined system the galaxy can be set up in a self- 
consistent way. By construction these models are in virial equi- 



librium. IKD95I have also demonstrated that their models are 
stable over many rotation periods. 

We chose model A of lKDQSl as a reference model (Tab.Q. 
The scale length of the disk isRi - 4 kpc. The total mass of the 
galaxy is 2.47 ■ lO" M0. The mass and radius of the disk are set 
to 0.34 • 10" Mq and 18.0 kpc. The corresponding parameters 
for the bulge are 0. 1 7 • 1 " Mq and 4.0 kpc and for the DM halo 
1 .94 -10" Mq and 87.9 kpc, respectively. The particle numbers 
for disk and bulge, A^d = 20000 and A^b = 10 000, are chosen 
to have particles of nearly equal mass. The DM halo is added 
as a static potential. The disk is initially in rotational equilib- 
rium with a small fraction of disordered motion superimposed. 
This results in a Toomre parameter Q = kctr/ (3. 3601.) exc eed- 
ing 1 .9 all over the disk. Thus, according to Polvachenko e t alJ 
ifl997D the disk is stable against all axisymmetric and non- 
axisymmetric perturbations. 

A cloudy gas phase is added by treating every fifth particle 
from the stellar disk as a cloud. The total mass in clouds is 
Meld ~ 6.9 • 10^ Mq and each cloud is randomly assigned a time 
of inactivity Ty^ between and 200 Myr. Initially, all clouds 
have the same mass (~ 1.7 ■ 10^ Mq) but a mass spectrum is 
built-up within the first Gyr of evolution. 

Finally, a diffuse gas component with a total mass of 
Mgas - 1 ■ 10^ Mq is added. This is more difficult because 
no simple equilibrium state exists for this gas phase: it is not 
only affected by the gravitational potential but also by cool- 
ing and heating processes. The latter varies locally and tempo- 
rally so that the diffuse gas can only reach a dynamical equi- 
libriu m as is also revealed in high- resolution simulations of the 
ISM jde Avillez & Breitschwerdtl [2004). Testing different ini- 
tial setups, starting far from equilibrium (for the diffuse com- 
ponent only), turned out to be most efficient. Then, the diffuse 
gas phase relaxes on a few dynamical timescales to a quasi- 
equilibrium state. In detail, we started with a slowly rotating 
homogeneous gaseous halo with a radius of 25 kpc and a tem- 
perature of 5.4 ■ 10"* K which relaxed to a quasi-equilibrium 
within a few 100 Myr 

Since we do not follow the chemical evolution of the 
galaxy, we apply solar abundances where a metallicity is 
needed. In order to treat the feedback of 'old' star particles 
each star particle is assigned an age. For the bulge stars the 
age is 10 Gyr and for disk stars - 10 Gyr. 

4.2. Evolution of the Reference Model 

Starting from these initial conditions the galaxy model is 
evolved for 3 Gyr with all processes switched on. After an ini- 
tial phase (500 Myr) dominated by the collapse of the gaseous 
halo the galaxy reaches a quasi-equilibrium defined by the in- 
terplay of the gravitational potential, radiative cooling, heat- 
ing by SNII and condensation and evaporation. When the equi- 
librium is reached the diffuse gas phase has formed a warm 
disk and a hot halo component. Their mass ratio of about 7 re- 
mains nearly constant for the rest of the simulation. The equi- 
librium state reached is almost independent of intial setup of 
the gas halo, i.e. the following evolution is not influenced by 
our choice. 
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Fig. 3. The evolution of stellar surface density (at 1 Gyr and 
3 Gyr). The stellar disk becomes slightly thicker with time but 
is otherwise stable with weak transient spiral patterns. The 
surface density plots (Fig. |3] and were compute d using the 
NEMO Stellar Dynamics Toolbox llTeubeiJll995l) . The parti- 
cles were binned in a grid with a cell size of 40 pc. The re- 
sulting map was then smoothed with a Gaussian kernel with 
FWHM = 1 kpc. 
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Fig. 4. The evolution of cloud surface density (at 1 Gyr and 
3 Gyr). The cloud disk is stable but the overall surface density 
is slowly decreasing. For details on how the surface density was 
computed see Fig.|3 



Snapshots of the evolution of the Milky Way model can be 
seen in Fig.|3lto|6l The stellar surface density plots in Fig.|3l 
demonstrate the stability of the disk during the 3 Gyr evolu- 
tion. Profiles of the surface density show that the radial mass 
distribution does not change except for the central 2 kpc where 
the surface density is increased by up to a factor of 5 caused 
by the dissipation-induced infall of clouds which later form 
stars. Face-on views show weak transient spiral patterns, e.g. 
at f = 2 Gyr. Edge-on views reveal a thickening of the disk. 
Throughout the disk the 95%-mass height zci5,sti increases by 
47% from 0.75 kpc to about 1.1 kpc. Only in the central 2 kpc 
can no thickening be seen. 

The evolution of the cloudy phase (Fig.|4} is influenced by 
dissipation and SF. The depletion of the cloudy medium by SF 
decreases its mean surface density. Due to inelastic cloud-cloud 
collisions the clouds also show a less smooth distribution. As a 
consequence a mass flow towards the center increases its cen- 
tral surface density. The radial profiles reveal that the surface 
density has decreased mostly around 5 kpc while it it almost 
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Fig. 5. The of evolution of gas density (at 1 Gyr and 3 Gyr). 
In each panel the xy- (bottom), xz -projection (top) and density 
scale are shown. Starting from a homogenous distribution the 
diffuse gas collapses and forms a thin gas disk within 1 Gyr. 
While the mass of the gaseous disk is nearly constant its ra- 
dius is growing slowly. Densities in the disk are of the order 
of 1 cm"^. The halo has densities from 10"^ cm"^ to 10"^ cm"^. 
Densities were computed on a lOOxlOO-grid using the SPH for- 
malism. 
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Fig. 6. The evolution of gas temperature (at 1 Gyr and 3 Gyr). 
In each panel the xy- (bottom), xz-projection (top) and temper- 
ature scale are shown. After the initial collapse the diffuse gas 
forms a disk at about 10^ K. The central part of the disk as well 
as small bubbles in the outer parts are heated by SNII to 10^ K. 
The halo gas it at 10*' - 10^ K. 

unchanged beyond 10 kpc. The cloud disk shows a thickening 
only outside of 1 1 kpc where the increase of zgs cid is very sim- 
ilar to that of the stellar disk. Inside that radius no thickening 
can be seen. In fact, due to dissipation z^s^m is even reduced 
inside 5 kpc by a factor of 2-3. 

The diffuse gas (Fig.|5j starts from a homogeneous distri- 
bution which collapses in less than 300 Myr to form a disk. 
Initially, the gaseous disk is only a few kpc in diameter but it 
is slowly growing in size. After 1 Gyr it has reached a diameter 
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Fig. 7. The evolution of the SFR with time (full line). The SFR 
decreases with time due to consumption of cloud mass. The av- 
erage SFR (dashed line) is 1 .6 Mq yr ' . The dotted line is the 
result of a simple analytic model accounting for the depletion 
of the cloudy medium by SF. 
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Fig. 10. The cloud mass spectrum. The different lines corre- 
spond to the mass spectrum after 1 Gyr (dotted), 2 Gyr (dashed) 
and 3 Gyr (full) of evolution. The fitted power-law mass func- 
tion is with a * -2.0 comparable to the observed mass spec- 
trum. 
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Fig. 8. The evolution of SF efficiency . SFE is up to 0.4 in the 
center of the galaxy. The average over the whole disk is about 
0.06. 
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Fig. 9. The open diamonds show the SFR per pc^ averaged 
over 100 Myr as a function of cloud surface densit y. The full 
line is a recent determination of the Schmidt Law by'KennicutJ 
(1998) with a slope of -1.4. The vertical dashed line indicates a 
cut-off density for SF in disk galaxies also found bv lKennicutll 
lll998l) . 



of 20 kpc. At 2 Gyr the diameter is ~ 30 kpc which is com- 
parable to that of the stellar disk. At the same time the mass 
of the gaseous disk is nearly constant. It has particle densities 
of 0.1 - 1 cm"^ and temperatures (Fig.|6} of less than lO'^K. 
Towards the center the gas reaches higher temperatures of up to 
several 10^ K. Some diffuse gas remains in the halo at densities 
of ~ 10 cm""* and less while its temperature is in the range of 
10^ - 10^ K. During the initial collapse the diffuse gas gains a 
significant amount of angular momentum transferred from the 
clouds by the drag force. For this reason the gaseous disk ro- 
tates with a rotation velocity of ~ 200 km s"'. Transient spiral 
patterns are also observed in the gas disk. Stellar feedback in- 
duces the emergence of bubbles of hot gas. These bubbles have 
diameters of up to 1 kpc and they cool within a few 10 Myr. 
The volume filling factor of the hot gas (r > 10"* K) in the disk 
is 30-40%. 

The velocity dispersion of the stars increases with time. At 
8 kpc the heating rate is about 1.3 kms ' Gyr"' and it is nearly 
constant over the whole disk. In order to see how much heat- 
ing is caused by unphysical two-body relaxation (which stems 
from the inability to resolve individual stars numerically) we 
compared the result with a pure N-body simulation and found 
that the heating rate in our full model is about 30% higher 
than in the pure N-body calculation. The velocity dispersion 
of the clouds decreases at 0.6 kms ' Gyr"' due to the energy 
dissipation by cloud-cloud collisions. As a result younger stars 
are born with a smaller velocity dispersion. This means that a 
part of the observed 'heating' in the age- velo city-relation might 
be ind eed due to cooling as suggested by iRocha-Pinto et alJ 

The average SFR is ~ 1 .6 Mq yr"' (Fig.0. Over the 3 Gyr 
of evolution the SFR decreases slowly from ~ 4.0MQyr"' to 
~ O.5M0yr"'. A simple analytic model indicates that this is 
due to the consumption of the clouds. The mean SFE is roughly 
constant at a level of e = 0.06. The maximum values of the lo- 
cal SFE are highest in the center with e = 0.2 - 0.4 (Fig. |8}. 
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It can be seen in Fig.|5]that the SFR also compares well to a 
recent determination of the Schmidt law by Kennicutt ( 1998). 
Kennicutt has also shown that the SF law in disk galaxies steep- 
ens abruptly below a threshold gas density of ~ IOMqpc"^. 
We also find in our simulations that the SF law steepens for 
lower gas densities. This is in agreement with Id CD models 
of the settling of a gas disk which show a similar threshold for 
the formation of a thin disk ( Burkert et al. 1992). Excluding the 
data below 10 M0 pc"^, we obtain a least-square fit of our data 
to Eq. (|6jl with w 1.7 + . 1 whi ch is in fair agreement with 
« * 1.4 + 0.15 of lKennicuttl ( ll998l) . 

The mass exchange by condensation and evaporation 
reaches a quasi-equilibrium after ~ 200 Myr. Before, evapo- 
ration is the dominant process due to the heating of the dif- 
fuse gas induced by the initial collapse of this component. In 
the later evolution, the condensation rate is decreasing from 
~ 2MQyr"' to O.lMQyr"'. The evaporation rate is slightly 
lower after 1 .5 Gyr due to the decreasing SN rate, resulting in 
a reduced heating of the diffuse ISM. 

The total cloud mass is slowly decreasing, showing the de- 
pletion by SF. The total mass of the diffuse gas remains al- 
most constant. A detailed look shows that more than half of 
the diffuse gas (65%) is located in the disk (R < 20 kpc and 
z < 1 kpc), while the rest is in the halo. 

In the initial model all clouds have the same mass. 
However, due to fragmentation and coalescence a stationary 
cloud mass spectrum is built-up within 1 Gyr (Fig. I10> . This 
mass spectrum can be fitted by a power-law with index a ^ -2 
with an uncertainty of 0.2. Once the mass spectrum is estab- 
lished it remains unchanged until the end of the simulation. 

5. Discussion 

In the previous section we presented a model for a Milky- Way- 
type galaxy. Since our description depends on several phase 
transitions and exchange processes we want to discuss their in- 
fluence in this section. In detail we discuss the SF, cloud mass 
spectrum and metallicities. 

The model has a number of parameters which can be cho- 
sen within the constraints of related observations or theory. 
We therefore performed simulations to investigate the influ- 
ence of some of those parameters using the model presented 
in Sect. 14.21 as a reference. We find that most models do not 
deviate strongly from the reference model. 

5.1. Star Formation 

The mean SFR is ~ 1.6Moyr"' in most models with a mean 
SF efficiency of ~ 6%. However, the SFR decreases in all mod- 
els from a few to below one solar mass per year. Observations 
on the other hand reveal a much more constant SF history 
jRocha-Pinto et aLi2000) . The decrease of the SFR can be re- 
produced in an analytic model taking into account the deple- 
tion of the cloudy medium by SF. An explanation for the dis- 
crepancy might be that our simulation does not include any in- 
fall to sustain a constant SFR. An infall rate of the order of 
~ 1 M0 yr"' would be sufficient for this provided the infalling 
gas can be used in the SF process. The necessary infall rate 



could be provid ed by high velocity clouds JBlitz etalJll99i 
Wakk er etalE9 99). 

The SF law is in good agreement with observations. Our 
data can be fitted to Eq. (|6|l with n ^ 1. 7 which is well 
within the observed range of n ^ 1.4- 2.0 llKennicuttlll9'95 
IWong & Blit3l2002t lioissier et alJl2003l) . It should be noted 
that this result is widely independent of the choice of our model 
parameters. Only models with a co nstant SF effic i ency s how a 
different behaviour with n - 1.2. iKoppen et alJ ( 1 19951) have 
shown that stellar feedback leads to a self-regulated SF follow- 
ing a Schmidt-type law. Our simulations confirm this result. 
We conclude that the observed index n of the SF law is a con- 
sequence of the self -regulation process which is enabled by al- 
lowing for a variable SEE. 

5.2. Cloud Mass Function 

The cloud mass function in our reference model has a slope 
of ff = -2. In most observations the mass spectrum is shal- 
lower with a ^ -1.5 (e.g. Sanders et al. 1985) though there 
are a few observations wit h larger values, e.g. a w -1.8 (e.g. 
Brand & Wouterloot 1995). Taking into account the uncertain- 
ties in both our model and observations we are able to resolve 
the cloudy phase of the ISM, i.e. each cloud particle in our sim- 
ulation represents an individual molecular cloud. 

Th eoretical models of coagulation by iField & Saslawl 
d 19651) also predict a - -1.5. These models include a frag- 
mentation by SF but (different to our model) SF occurs only 
when clouds reach a maximum mass. The resulting fragments 
are then added to the lowest mass bin. In a more similar model, 
which also takes into account a mass-radius relation (Eq. 
lElme greerii C 1 9 8 9 ) found a mass spectrum with a - -2, show- 
ing that our simulations can reproduce these results. 

5.3. Metallicity 

We compared two models with metallicities Z = 1/2 Zq and 
Z - 1.5 Zq. The change of metallicity is only moderate and 
it does not affect the dynamical evolution of the galaxy sig- 
nificantly. This result verifies that using a constant metallicity 
is a valid approximation on short time scales. Using a simple 
estimate we find that for 3 Gyr the expected change in metallic- 
ity would be of the order I/IOZ©, i.e. well within in the range 
given by the two test models. Therefore, a model including stel- 
lar yields would not change the results as long as the late evo- 
lution (i.e. the last few Gyr) of galaxies is simulated. However, 
stellar yields would have to be included when the early forma- 
tion stages are considered. 

5.4. Tine Muiti-Phase ISM 

The ISM in our models can be divided into three phases: molec- 
ular clouds and diffuse components, which can be either warm 
or hot. The total gas mass fraction after 3 Gyr of evolution is 
about 7% in all our models. The ratio of the mass of clouds to 
the mass of diffuse gas is mainly determined by condensation 
and evaporation. If the evaporation rate is higher, for exam- 
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pie due to higher feedback energy this ratio is closer to unity. 
80 - 90% of the diffuse gas is found in the warm gaseous disk. 
Only ~ 15% remain in the halo in a hot gas phase. In the low- 
metallicity model this fraction is larger by a factor of 2 due 
to less efficient cooling. Samland & Gerhard (2003) have also 
presented a model of a disk galaxy. In their model the ratio of 
cloud and diffuse gas mass is 3 which is comparable to our re- 
sults. They also find a hot g aseous halo. In a similar calculation 
LSemelin & Co mbes' ( 2002) find equal masses in clouds and gas 
independent of the feedback energy, but their models do not in- 
clude phase transitions between the gas phases by condensation 
and evaporation. 

We have found a filling facto r f of 0.3 - 0.4 for the hot gas 
in the disk. McKee & Ostriker ( 19771) predicted 0.7, a larger 
value, but ou r value is in better agreement with the more recent 
findings of Norm an & Ikeuchil fl989i) who predicted a filling 
factor of about 0.2. Furthermore, th ree dimensional simulation 
of a region of the galactic disk bv Ide AvillezI ll200 0) show a 
filling factor for gas with T > 10"* K of ~ 0.4 in good agreement 
with our results. We find an increase of 0.2 in / if the energy 
injection by S NII is higher by a factor 4 which is comparable 
to the results ofl de Avillez & Breitschwerdtl (l2004h . 

One shortcoming in our treatment of the ISM is the lack of 
a detailed cooling below 10"^ K. As a consequence we do not al- 
low clouds to form directly by thermal instabilities. However, 
due to low mean densities and the turbulent state of the dif- 
fuse gas we expect the formation of clouds in expanding SN 
shells (which is included) to be much more important. Again, a 
more sophisticated treatment of thermal instabilities would be 
required for modelling the early evolution of (gas rich) galax- 
ies. 

5.5. Live DM Halo 

So far, in all presented models the DM halo was treated as a 
static potential. For the case of the late evolution of an isolated 
galaxy a static halo is reasonable. However, DM particles have 
to be used, e.g., for interacting galaxies or early galactic evo- 
lution. In one model we therefore used a DM halo made of 
particles in order to test how a live halo changes our results. 
The DM halo was modelled with 10^ particles, thus DM parti- 
cles are twice as massive as the initial masses of star particles. 
As expected we do not find significant changes in the evolution 
when using a live DM halo. The CPU-time needed in this run 
is higher by a factor of 2 (lOcpu-days instead of 5). 

6. Summary and Conclusions 

We have presented a new particle-based code for modelling 
the evolution of galaxies. In addition to the collisionless stellar 
and DM components it has incorporated a multi-phase descrip- 
tion of the ISM combining standard SPH with a sticky particle 
method to describe diffuse warm/hot gas and cold molecular 
clouds respectively. We have included a new description for SF 
with a local SF efficiency depending on the local pressure of the 
ISM and the mass of star forming clouds. Stellar life times and 
an IMF are taken into account as well as SNII and PNe. The 



gas phases can mix by means of condensation and evaporation 
and clouds are subject to a drag force. 

We used this description to model 3 Gyr of evolution of a 
Milky-Way-type galaxy. Overall, the evolution of our model 
galaxy is quite stable. Furthermore, a parameter study showed 
that the evolution does not change significantly when varying 
model parameters like feedback energy, cloud radius, SF time 
scale or metallicity in a reasonable range. 

Our models produce a SF law with an index n - 1.7 in 
good agreement with observations. For models with a constant 
SF efficiency e the index is with n = 1 .2 smalle r than obser- 
vations suggest. We conclude, in agreement with Kop pen et alJ 
( 1995), that the SF law is the result of a self-regulation process 
enabled by our new treatment of SF. The average SFR is also 
comparable to observations. Over time the SFR decreases by 
a factor of 4-8 because we have not included any gas infall. 
To get a more constant SFR a gas infall should be taken into 
account at a rate of ~ 1 Mq yr"' . The cloud mass spectrum in 
our models (with a = -2) is a little bit steeper than most ob- 
servation suggest. Nonetheless, each cloud particle in our code 
essentially represents an individual molecular cloud. 

In all models most (85%) of the diffuse gas forms a warm 
(~ 10"* K) rotating disk. A hot gaseous halo with temperatures 
of lO*"- 10^ K and densities of 10"^ 

cm has been built up. Hot 
gas is also found in the disk in bubbles accompanied by SNII. 
The filling factor of the hot disk gas is about 0.4 comparing 
well to observations as well as models of the ISM. This result 
could be used to find more realistic conditions for the initial 
model. 

In a next step we will apply our model to interacting galax- 
ies in order to follow the evolution of the ISM and the SF 
history in perturbed galactic systems. The chemical evolution 
should be embedded in order to allow for a more detailed con- 
frontation with observations. 
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Appendix A: Determination of the star formation 
efficiency 



lEE97l derived an implicit relation between cloud mass niM, 
pressure Pg^^and SEE e (M5 = 10^ Mq) 



and 



Wield 

M5 



3 . /'Weld 



-1/4 /p 

1/4 



-5/8 



-5/8 



1 + 



(6/^)2/5-1.4 



0.56 



(A.2) 



4^ , when- > 1 



lEE97l use Pq - 3 ■ 10 cm K for the pressure in the solar 
neighbourhood. ^ is defined as 



. /'«eld\ /^^gas\ 



(A.3) 



The constants Aq = 180 and ^0 - 0.1 are determined from 
observations of star forming clouds in the solar neighbourhood. 

In each SF event Eq. ( lA.U and Eq. ( IA.2t are solved nu- 
merically to determine e. The resulting SEE is highest for high 
mass clouds in a high pressure environment. Typical SEE are 
between 5 and 10%. 

Appendix B: SPI-I implementation 

The idea of SPH is that a fluid is represented by particles. 
Any physical property A of the fluid at the position r^, is evalu- 
ated by means of a locally weighted average: 



Pi 



-W(r,,„Ji). 



(B.l) 



The summation is performed over particles q with mass rriq, 
density pq. Aq is the value of A at and the distance Tpq is 
Tq. We apply th e spherically symmetric smooth ing kernel 



W(r, h) suggested bv lMonaghan & Lattanziol ( 1 19851 



W(r, h) = 



n9 



forO<f 

2\hl 4\h) * 
4\ h) ' 



< 1 



I 0, 



for 1 < f < 2 



for ^ > 2. 



(B.2) 
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A variable smoothing length h keeping the number A^sm of near- 
est neighbours nearly constant at about 32 for all SPH particles 
is applied. 

The equation of motion of an SPH particle p is given by 



dfi 



= - — VPp+af'-V^p, 
Pp 



(B.3) 



where O,, is the gravitational potential, Pp the pressure and a™*^ 
an artificial viscosity term to handle shocks. We use a symmet- 
ric expression to estimate (VP)/p: 



Hgas and can be easily rescaled. Starting from a given temper- 
ature the new temperature after a (rescaled) time step is then 
simply looked up from the table. This provides a fast way of 
computing the change in temperature, i.e. the cooling rate, for 
any given temperature, density and time step. The cooling rate 
A for each particle is updated each individual time step. 

A more detailed description of the SPH method can be 
found in .Monaghan ( 1992) . 



(B.4) 



where Wpq indicates that the arithmetic mean of the two 
kernel function for the individual smoothing lengths hp and 
hn is used. For the artific ial viscosity term a]^"*'^ we use 
llningold&Monaghanll983l) 

af'^-Y,mqnp,VWp„ (B.5) 



with 



n 



pq 



Ppq 



(B.6) 



where ppg and Cpq are the arithmetic means of density and 
sound speed of particles p and q. fipq is defined by 



t^pq = 



hpq{rlq/hlg + r]l) 
forv 



for \pq -Tpg <0 



pq Tpq > 



(B.7) 



where ypq is the relative velocity of particle p and q. The vis- 
cosity parameters are set to - 2, [5^ - I and tj^ - 0.1. The 
evolution of the specific internal energy u is computed by 



dup 
dt 



Pi 



+ n 



pq 



'pq 



r- A 



(B.8) 



r and A account for heating and cooling processes. 

The time integration is done by a leap-frog integrator. 
Individual time steps are used for each particle. The time step 
for an SPH particle p is chosen to satisfy a modified form of 
the Courant condition as suggested b v Monaghan (. 1992) 



Mp < C 



hp |V ■ Vp| +Cp + 1.2 [uyCp +/3y maxg |jUp,^|) 



(B.9) 



where C = 0.4 is the Courant number, hp the smoothing length, 
Cp the sound speed, \p the velocity of particle p. 

Because cooling time scales are usually much shorter than 
the dynamical times, impractically short time steps might 
occur. Therefore, Eq. ( IB.8> is often solved semi-implicitly. 
However, we chose a different approach: a fourth-order Runge- 
Kutta scheme is used to integrate the cooling function before- 
hand with a sufficiently small time step and a look-up table with 
times and temperatures is created. The number density Wgas is 
fixed to one value for this but cooling times scale linear with 



